% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%     case1： 幂函数f(x) = sqrt(X)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% X = 1:100;
% Y = sqrt(X);
% % 可视化
% figure(1)
% plot(X, Y, 'y');
% title('函数sqrt(X)图像');
% xlabel('x');
% ylabel('sqrt(x)');
% % saveas(1, '函数sqrt(X)图像.jpg');
% 
% t = X;
% h = 1;
% % 使用求导公式求导
% val_1 = 0.5 * (1 ./ sqrt(X));
% % 使用插值型函数求导
% val_2 = Inter_diff(X, Y, t, h);
% % 两种方式的差值
% D_val = abs(val_1 - val_2);
% % D_val
% % 可视化
% figure(2)
% plot(X, val_1, 'r', X, val_2, 'g');
% title('函数sqrt(X)使用公式法、本算法求导的对比')
% hold on 
% stem(X, D_val, 'ob');
% legend('公式法','本算法', '相对误差');
% % saveas(2, '函数sqrt(X)对比.jpg');
% hold off
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%     case2：多项式函数f(x) = x^7 + x^4 + x^2 + 3x +1
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% X = 1:100;
% Y = X.^7 + X.^4 + X.^2 + 3 .* X + 1;
% % 可视化
% figure(3)
% plot(X, Y, 'y');
% title('函数x^7+x^4+x^2+3x+1图像');
% xlabel('x');
% ylabel('x^7+x^4+x^2+3x+1');
% % saveas(3, '函数x^7+x^4+x^2+3x+1图像.jpg');
% 
% t = X;
% h = 1;
% % 使用求导公式求导
% val_1 = 7 .* X.^6 + 4 .* X.^3 + 2 .* X + 3;
% % 使用插值型函数求导
% val_2 = Inter_diff(X, Y, t, h);
% % 两种方式的差值
% D_val = abs(val_1 - val_2);
% % D_val
% % 可视化
% figure(4)
% plot(X, val_1, 'r', X, val_2, 'g');
% title('函数x^7+x^4+x^2+3x+1使用公式法、本算法求导的对比')
% hold on 
% stem(X, D_val, 'ob');
% legend('公式法','本算法', '相对误差');
% % saveas(4, '函数x^7+x^4+x^2+3x+1对比.jpg');
% hold off
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%     case3：三角函数f(x) = sin(x)+cos(2x)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% X = linspace(-3 * pi, 3 * pi, 100);
% Y = sin(X) + cos(2 .* X);
% % 可视化
% figure(5)
% plot(X, Y, 'y');
% title('函数sin(x)+cos(2x)图像');
% xlabel('x');
% ylabel('sin(x)+cos(2x)');
% % saveas(5, '函数sin(x)+cos(2x)图像.jpg');
% 
% t = X;
% h = (6 * pi) / 100;
% % 使用求导公式求导
% val_1 = cos(X) - 2 * sin(2 .* X);
% % 使用插值型函数求导
% val_2 = Inter_diff(X, Y, t, h);
% % 两种方式的差值
% D_val = abs(val_1 - val_2);
% % D_val
% % 可视化
% figure(6)
% plot(X, val_1, 'r', X, val_2, 'g');
% title('函数sin(x)+cos(2x)使用公式法、本算法求导的对比')
% hold on 
% stem(X, D_val, 'ob');
% legend('公式法','本算法', '相对误差');
% % saveas(6, '函数sin(x)+cos(2x)对比.jpg');
% hold off
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%     case4：分段f(x) = x^3，x < 0; x^2，x >= 0
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% X = -99:1:100;
% Y = zeros(1,200);
% for i = 1:200
%     if X(i) < 0
%         Y(i) = X(i)^3;
%     else
%         Y(i) = X(i)^2;
%     end
% end
% 
% % 可视化
% figure(7)
% plot(X, Y, 'y');
% title('分段函数图像');
% % saveas(7, '分段函数图像.jpg');
% 
% t = X;
% h = 1;
% % 使用求导公式求导
% val_1 = zeros(1, 200);
% for i = 1:200
%     if X(i) < 0
%         val_1(i) = 3 * X(i)^2;
%     else
%         val_1(i) = 2 * X(i);
%     end
% end
% % 使用插值型函数求导
% val_2 = Inter_diff(X, Y, t, h);
% % 两种方式的差值
% D_val = abs(val_1 - val_2);
% D_val
% % 可视化
% figure(8)
% plot(X, val_1, 'r', X, val_2, 'g');
% title('分段函数使用公式法、本算法求导的对比')
% hold on 
% stem(X, D_val, 'ob');
% legend('公式法','本算法', '相对误差');
% % saveas(8, '分段函数对比.jpg');
% hold off

syms x
y = (exp(x^2) - 8 * x^4) / sin(x);
y_1 = diff(y,x,1);
res = double(subs(y_1, x, 0.75));

X = 0.5:0.05:1.5;
Y = (exp(X.^2) - 8 * X.^4) ./ sin(X);
t = 0.75;
h = 0.05;
val = Inter_diff(X, Y, t, h);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 函数功能：  插值求导
% 参数意义：  1.X：插值点的横坐标数组
%            2.Y：插值点的纵坐标数组
%            3.t：需要求导数的点数组
%            4.h：插值点的步长
% 返回值：    value——待求导点数组的导数值
% 注意：函数只接受等间距插值的数据
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function value = Inter_diff(X, Y, t, h)
% 判断输入数据的维数是否符合
if (length(X) ~= length(Y))
    error('The input data dimensions are different! Please check the input parameters.');
elseif ((isempty(X) == 1) || (length(X) == 1))
    error('The input data dimension must be at least 2 !');     
end

n = length(X);
k = length(t);
% 构造导数值矩阵
value = zeros(1, k);
% 根据参数X长度选择不同的插值求导方式
switch n
    case 2
        for i = 1:k
            % 选择两点公式
            value(i) = point2(X, Y);
        end
    case 3
        for i = 1:k
            % 选择三点公式
            value(i) = point3(X, Y, h, t(i));
        end
    case 4
        % 由于参数X长度大于3小于5，只能选择三点公式
        for i = 1:k
            index = 0;
            % 获取带求求导点在X中的索引
            index = find(X == t(i));
            % 如果索引不是首尾元素，选择左右邻点构成三个点，再调用三点公式
            % 如果是首元素，选择前三个点；如果是尾元素，选择最后三个点
            if (((index - 1) >= 1) && ((index + 1) <= n))
                % 选择插值点
                X_new = X(1, index-1:index+1);
                Y_new = Y(1, index-1:index+1);
                % 选择三点公式
                value(i) = point3(X_new, Y_new, h, t(i));
            elseif ((index - 1) < 1)
                % 选择插值点
                X_new = X(1, index:index+2);
                Y_new = Y(1, index:index+2);
                % 选择三点公式
                value(i) = point3(X_new, Y_new, h, t(i));
            else
                % 选择插值点
                X_new = X(1, index-2:index);
                Y_new = Y(1, index-2:index);
                % 选择三点公式
                value(i) = point3(X_new, Y_new, h, t(i));
            end
        end
    case 5
        for i = 1:k
           % 选择五点公式
           value(i) = point5(X, Y, h, t(i)); 
        end
    otherwise
        % 由于参数X长度大于5，直接选择五点公式求导
        for i = 1:k
           index = 0;
           % 获取带求求导点在X中的索引
           index = find(X == t(i));
           % 如果索引不是最前面两个和最后面两个元素，左右邻点各选择两个构成五个点，再调用三点公式
           % 如果是最前面两个元素，选择前五个点；如果是最后两个元素，选择最后五个点
           if (((index - 2) >= 1) && ((index + 2) <= n))
               % 选择插值点
               X_new = X(1, index-2:index+2);
               Y_new = Y(1, index-2:index+2);
               % 选择五点公式
               value(i) = point5(X_new, Y_new, h, t(i));
           elseif ((index - 2) < 1)
               % 选择插值点
               X_new = X(1, 1:5);
               Y_new = Y(1, 1:5);
               % 选择五点公式
               value(i) = point5(X_new, Y_new, h, t(i));
           else
               % 选择插值点
               X_new = X(1, n-4:n);
               Y_new = Y(1, n-4:n);
               % 选择五点公式
               value(i) = point5(X_new, Y_new, h, t(i));
           end
        end
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 函数功能：  两点求导
% 参数意义：  1.X：插值点的横坐标数组
%            2.Y：插值点的纵坐标数组
% 返回值：    value——待求导点数组的导数值
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function value = point2(X, Y)
% 运用两点公式
value = (Y(2) - Y(1)) / (X(2) - X(1));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 函数功能：  三点求导
% 参数意义：  1.X：插值点的横坐标数组
%            2.Y：插值点的纵坐标数组
%            3.t：需要求导数的点数组
%            4.h：插值点的步长
% 返回值：    value——待求导点数组的导数值
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function value = point3(X, Y, h, t)
% 根据三点公式分别求取各点的导数值
% f'(x1)、f'(x2)、f'(x3)计算如下
f_1 = (-3 * Y(1) + 4 * Y(2) - Y(3)) / (2 * h);
f_2 = (Y(3) - Y(1)) / (2 * h);
f_3 = (Y(1) - 4 * Y(2) + 3 * Y(3)) / (2 * h);
% 构建导数矩阵
f = [f_1 f_2 f_3];
% 获取待求点t在X中的索引
index = find(X == t);
% 根据indx从f找到导数值
value = f(index(1));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 函数功能：  五点求导
% 参数意义：  1.X：插值点的横坐标数组
%            2.Y：插值点的纵坐标数组
%            3.t：需要求导数的点数组
%            4.h：插值点的步长
% 返回值：    value——待求导点数组的导数值
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function value = point5(X, Y, h, t)
% 根据五点公式分别求取各点的导数值
% f'(x1)、f'(x2)、f'(x3)、f'(x4)、f'(x5)计算如下
f_1 = (-25 * Y(1) + 48 * Y(2) - 36 * Y(3) + 16 * Y(4) - 3 * Y(5)) / (12 * h);
f_2 = (-3 * Y(1) - 10 * Y(2) + 18 * Y(3) - 6 * Y(4) + Y(5)) / (12 * h);
f_3 = (Y(1) - 8 * Y(2) + 8 * Y(4) - Y(5)) / (12 * h);
f_4 = (-Y(1) + 6 * Y(2) - 18 * Y(3) + 10 * Y(4) + 3 * Y(5)) / (12 * h);
f_5 = (3 * Y(1) - 16 * Y(2) + 36 * Y(3) -48 * Y(4) + 25 * Y(5)) / (12 * h);
% 构建导数矩阵
f = [f_1 f_2 f_3 f_4 f_5];
% 获取待求点t在X中的索引
index = find(X == t);
% 根据indx从f找到导数值
value = f(index(1));
end